rm(list=ls())
library(arfima)
library(matrixStats)
library(lattice)
#
# Create Function
#

phi <- .2
theta <- .21
d <- .2
i <- 0
x <- 100
nsim <-10

experiment <- function(phi, theta, d, i, x,nsim)
{
  res <- vector(mode="list", length=nsim)
  for(m in 1:nsim){
    # Simulate
    ts <- arfima.sim(x,model = list(theta = theta,
                                    phit = phi,
                                    dfrac=d, 
                                    dint =i))
    # Fit and Summarize
    arfm <- summary(arfima(ts,order = c(1, 1, 1),dmean=FALSE))
    arfm2 <- arfm[5]
    arfm3 <- arfm2$coef[[1]]
    
    #
    res[[m]] <- c(mod = arfm2,d)
  }
  class(res) <- "arfsim"
  return(res)
}


#
# Declare values
#

phi <- c(.5)
theta <- c(.3)
d <- c(-.2)
i <- c(1)
x <- c(50,100,250,500,1000,1500)
nsim <- 100

eg <- expand.grid(x,d,phi,theta,i)
eg  
colnames(eg) <-c("x","d","phi","theta","i")
eg

allexp <- vector(mode="list", length=nrow(eg))

#
# Run Experiments
#


for(v in 1:nrow(eg)){
  allexp[[v]] <- experiment(x =eg[v,1],
                            d=eg[v,2],
                            phi=eg[v,3],
                            theta=eg[v,4],
                            i=eg[v,5],
                            nsim=nsim)
}








#
# Summarize Results
#

# Output in a matrix (eg,sim)

output<- as.matrix(allexp)
output.2 <- do.call(rbind,allexp)
output.3 <- apply(output.2,c(1,2),function(x)x[[1]]$mod.coef[1])
output.4 <- matrix(lapply(output.3,function(x)x[[1]]),nrow=nrow(eg),ncol=nsim)

output.4[[2]]

matrix.index <- expand.grid(x,d)
matrix.index

# Matrix of D estimates
output.5 <- matrix(as.numeric(lapply(output.4,function(x)x[3])),nrow=nrow(eg),ncol=nsim)
rownames(output.5) 
colnames(output.5) <- (1:nsim)

# Matrix of D Standard Errors
output.6 <- matrix(as.numeric(lapply(output.4,function(x)x[6])),nrow=nrow(eg),ncol=nsim)
rownames(output.6) 
colnames(output.6) <- (1:nsim)

# Matrix of D
output.d <- apply(output.2,c(1,2),function(x)x[[1]][[2]])
rownames(output.d) 
colnames(output.d) <- (1:nsim)



#################
## D Estimates ##
#################

#d=.0#

d.50.8 <- output.5[1,]                  
d.100.8 <- output.5[2,]  
d.250.8 <- output.5[3,]
d.500.8 <- output.5[4,]
d.1000.8 <- output.5[5,]
d.1500.8 <- output.5[6,]

e.50.8 <- output.6[1,]                  
e.100.8 <- output.6[2,]  
e.250.8 <- output.6[3,]
e.500.8 <- output.6[4,]
e.1000.8 <- output.6[5,]
e.1500.8 <- output.6[6,]

#Average Estimates
rowMeans(output.5)

#################
## Average CIs ##
#################

#
# d.8
#

#50

d.low.50.8 <- c(d.50.8 - 2.009*(e.50.8))
d.high.50.8 <- c(d.50.8 + 2.009*(e.50.8))
d.is.50.8 <- cbind(d.low.50.8,d.50.8,d.high.50.8)
d.is.50.8
l.1.50.8 <- colMeans(d.is.50.8)

#100

d.low.100.8 <- c(d.100.8 - 1.984*(e.100.8))
d.high.100.8 <- c(d.100.8 + 1.984*(e.100.8))
d.is.100.8 <- cbind(d.low.100.8,d.100.8,d.high.100.8)
d.is.100.8
l.2.100.8 <- colMeans(d.is.100.8)

#250

d.low.250.8 <- c(d.250.8 - 1.96*(e.250.8))
d.high.250.8 <- c(d.250.8 + 1.96*(e.250.8))
d.is.250.8 <- cbind(d.low.250.8,d.250.8,d.high.250.8)
d.is.250.8
l.3.250.8 <- colMeans(d.is.250.8)

#500

d.low.500.8 <- c(d.500.8 - 1.96*(e.500.8))
d.high.500.8 <- c(d.500.8 + 1.96*(e.500.8))
d.is.500.8 <- cbind(d.low.500.8,d.500.8,d.high.500.8)
d.is.500.8
l.4.500.8 <- colMeans(d.is.500.8)


#1,000

d.low.1000.8 <- c(d.1000.8 - 1.96*(e.1000.8))
d.high.1000.8 <- c(d.1000.8 + 1.96*(e.1000.8))
d.is.1000.8 <- cbind(d.low.1000.8,d.1000.8,d.high.1000.8)
d.is.1000.8
l.5.1000.8 <- colMeans(d.is.1000.8)


#1,500

d.low.1500.8 <- c(d.1500.8 - 1.96*(e.1500.8))
d.high.1500.8 <- c(d.1500.8 + 1.96*(e.1500.8))
d.is.1500.8 <- cbind(d.low.1500.8,d.1500.8,d.high.1500.8)
d.is.1500.8
l.6.1500.8 <- colMeans(d.is.1500.8)


t.o.8 <- cbind(l.1.50.8,l.2.100.8,l.3.250.8,l.4.500.8,l.5.1000.8,l.6.1500.8)
table.out.8 <- t(t.o.8)
rownames(table.out.8) 
colnames(table.out.8) <- c("d.low","d.est","d.high")
table.out.8

# Plots of Distributions
par(mfrow=c(1,1))
plot (density(d.50.8), col="red",ylim=c(0,6),xlim=c(-1,1)) 
lines (density(d.100.8), col="green") 
lines (density(d.250.8), col="blue")
lines (density(d.500.8), col="orange")
lines (density(d.1000.8), col="black")
lines (density(d.1500.8), col="grey")
abline(v=-.2)
